/*
Output: Figure SI1. Estimated Coefficients Relating Owner and Parcel Variables to Percent of Area Harvested
*/

run "/Users/holyfantastic/Dropbox/Forest/Paper/Nature/FinalDocuments/Code/000declare_path" // change to your local path

set scheme white_jet 

use "ParcelLevelDataSept282023_VerNov82024_R1.dta" , clear

{

la var forestAcres "Forested acres of parcel"
la var forestAcresTotal "Total forested acres owned by owner"

*Requiring at least one acre of forest
drop if forestAcres<1
drop if taxidnum == ""

*Focusing on pre enrollment years (pre-2020)
drop if year>2019

*Mean harvesting by year
bysort year: egen pcut_mean=mean(Landtr_Cut_mean_forest)

*clean variables
gen forestAcres2=forestAcres^2
gen distance2=distance^2
gen long_owner = YearsOwned > 10


*recode ndvi
replace ndvi_landtr_mean_forest = . if ndvi_landtr_mean_forest <=0 
replace ndvi_landtr_mean_forest = ndvi_landtr_mean_forest / 1000

*demean ndvi and ndvi2

foreach var in ndvi_landtr_mean_forest  {
	qui:sum `var'
	replace `var' = `var' - r(mean)
}
gen  ndvi_landtr_mean_forest2 = ndvi_landtr_mean_forest^2
gen  ndvi_landtr_mean_forest3 = ndvi_landtr_mean_forest^3

// gen  ndvi_landtr_mean_forest2_dm = ndvi_landtr_mean_forest^2

*Let's do everything with the percent harvested, not the share harvested
replace Landtr_Cut_mean_forest = Landtr_Cut_mean_forest * 100

**# regression analysis

**## look at key variables

tab group
egen taxidnum_id  = group(taxidnum)
cap egen county_id = group(county)
xtset  taxidnum_id year
sum Landtr_Cut_mean_forest

}


**## Coefficient figure

global controls  i.forestHaGroup i.roadGroup  ///
				i.yearGroup i.iPerson  ///
				i.iLocal
				
reghdfe Landtr_Cut_mean_forest i.group l.ndvi_landtr_mean_forest l.ndvi_landtr_mean_forest2  $controls  if year >= YearBought	 , a(i.county_id##i.year) vce(cl taxidnum_id)



// preserve 

matrix A=J(20,3,.)
local k = 1
	
{
// 	/*Gap line*/
// 	mat A[`k',1]=`k'
// 	mat A[`k',2]=.
// 	mat A[`k',3]=.
// 	local k = `k' + 1
	
	*forestHaGroup 
	mat A[`k',1]=`k'
	mat A[`k',2]=_b[2.forestHaGroup]
	mat A[`k',3]=_se[2.forestHaGroup]
	local k = `k' + 1
	
	mat A[`k',1]=`k'
	mat A[`k',2]=_b[3.forestHaGroup]
	mat A[`k',3]=_se[3.forestHaGroup]
	local k = `k' + 1
	
	mat A[`k',1]=`k'
	mat A[`k',2]=_b[4.forestHaGroup]
	mat A[`k',3]=_se[4.forestHaGroup]
	local k = `k' + 1

	/*Gap line*/
	mat A[`k',1]=`k'
	mat A[`k',2]=.
	mat A[`k',3]=.
	local k = `k' + 1
	mat A[`k',1]=`k'
	mat A[`k',2]=.
	mat A[`k',3]=.
	local k = `k' + 1
	
	*roadGroup 
	mat A[`k',1]=`k'
	mat A[`k',2]=_b[2.roadGroup]
	mat A[`k',3]=_se[2.roadGroup]
	local k = `k' + 1

	mat A[`k',1]=`k'
	mat A[`k',2]=_b[3.roadGroup]
	mat A[`k',3]=_se[3.roadGroup]
	local k = `k' + 1
	
	mat A[`k',1]=`k'
	mat A[`k',2]=_b[4.roadGroup]
	mat A[`k',3]=_se[4.roadGroup]
	local k = `k' + 1
	
	/*Gap line*/
	mat A[`k',1]=`k'
	mat A[`k',2]=.
	mat A[`k',3]=.
	local k = `k' + 1
	mat A[`k',1]=`k'
	mat A[`k',2]=.
	mat A[`k',3]=.
	local k = `k' + 1
	
	*yearGroup 
	mat A[`k',1]=`k'
	mat A[`k',2]=_b[2.yearGroup]
	mat A[`k',3]=_se[2.yearGroup]
	local k = `k' + 1

	mat A[`k',1]=`k'
	mat A[`k',2]=_b[2.yearGroup]
	mat A[`k',3]=_se[2.yearGroup]
	local k = `k' + 1
	
	mat A[`k',1]=`k'
	mat A[`k',2]=_b[4.yearGroup]
	mat A[`k',3]=_se[4.yearGroup]
	local k = `k' + 1
	
	/*Gap line*/
	mat A[`k',1]=`k'
	mat A[`k',2]=.
	mat A[`k',3]=.
	local k = `k' + 1
	mat A[`k',1]=`k'
	mat A[`k',2]=.
	mat A[`k',3]=.
	local k = `k' + 1
	
	*iPerson
	mat A[`k',1]=`k'
	mat A[`k',2]=_b[1.iPerson]
	mat A[`k',3]=_se[1.iPerson]
	local k = `k' + 1	
	
	/*Gap line*/
	mat A[`k',1]=`k'
	mat A[`k',2]=.
	mat A[`k',3]=.
	local k = `k' + 1
	mat A[`k',1]=`k'
	mat A[`k',2]=.
	mat A[`k',3]=.
	local k = `k' + 1
	
	*iLocal
	mat A[`k',1]=`k'
	mat A[`k',2]=_b[1.iLocal]
	mat A[`k',3]=_se[1.iLocal]
	local k = `k' + 1
	
	/*Gap line*/
	mat A[`k',1]=`k'
	mat A[`k',2]=.
	mat A[`k',3]=.
	local k = `k' + 1
	
}
	

	mat list A
	clear
	svmat A
	sum A1 
	drop if A1>r(max)
	drop if A2 == .
	keep A1 A2 A3
	
	rename A2 coef
	rename A3 se
	rename A1 n1
	gen up = coef + 1.96 * se
	gen low = coef - 1.96 * se
	replace n1 = - n1
	replace n1 = n1 / 2	
	global tmin = -1
	global color "gs7"
	sum n1
	global ymin = r(min)
	sum low 
	global xmin = min(r(min),$tmin)
	sum up 
	global xmax = r(max)
	global step = round(($xmax - $xmin)/5 ,0.001)
	
	
	set scheme white_jet 
	
	
	
	**# plot
	* Figure: Heterogeneity 
	*------------------------------------------------------begin------------
	#delimit ;
	twoway 
		   (scatter n1 coef , msize(small)  mc($color) legend(off))
			(rcap up low n1 , lp(dash) lc($color) hori  legend(off) lw(medthick)) , 	
		  	 ytitle("")
/* 		   ylabel($ymin($ymin)0 0 " " $ymin " ") */
 		   ylabel(-10(10.5)0.5 0.5 " " -10 " ") 
		   xlabel($xmin(0.5)$xmax , grid) 
		   xline(0, lp(dash))
		   xtit("Effects of categorical variables on harvest, Coef. and 95% CIs" ,size(small)) 
		text(0 $xmin "Forest area (Ha)",place(e) size(small)) 
			text(-0.5 $xmin "          4-20 ha",place(e) size(vsmall) color($color)) 
			text(-1 $xmin "          40-79 ha",place(e) size(vsmall) color($color)) 
			text(-1.5 $xmin "          80 or more ha",place(e) size(vsmall) color($color)) 
		  text(-2.5 $xmin "Distance to road",place(e) size(small))
			text(-3 $xmin "          100-249 meters",place(e) size(vsmall) color($color)) 
			text(-3.5 $xmin "          250-500 acress",place(e) size(vsmall) color($color)) 
			text(-4 $xmin "          >500 meters",place(e) size(vsmall) color($color)) 
		  text(-5 $xmin "Years owned",place(e) size(small))
			text(-5.5 $xmin "          5-9 years",place(e) size(vsmall) color($color)) 
			text(-6 $xmin "          10-14 years",place(e) size(vsmall) color($color)) 
			text(-6.5 $xmin "          15 or more years",place(e) size(vsmall) color($color)) 
		  text(-7.5 $xmin "Owner is a person, not entity",place(e) size(small) )
			text(-8 $xmin "          Yes",place(e) size(vsmall) color($color))
		  text(-9 $xmin "Owner's address is within 40 km of parcel",place(e) size(small) )
			text(-9.5 $xmin "          Yes",place(e) size(vsmall) color($color))
			graphregion(margin(small))
	;
	#delimit cr
	*------------------------------------------------------over------------

foreach fig_path in figure {
	
	 dis  "$`fig_path'/fig_coef_plot_R1"
		
	graph export "$`fig_path'/fig_coef_plot_R1.png", as(png) name("Graph") replace
	graph export "$`fig_path'/fig_coef_plot_R1.pdf", as(pdf) name("Graph") replace

}
